%  This simulates the posterior for the Stock-Watson-Romer UC model with
%  stochastic volatility and measurement error.  This version is for US
%  wholesale price inflation.  A post WWII overlap between clean and noisy
%  measures is exploited to identify measurement error parameters.

% This version has AR(1) measurement error.

clear
NG = 100000; % number of draws from Gibbs sampler per data file
NF = 40;

% catalog of data files
DFILE(1,:) = ['swuc_wpi_AR1_01'];
DFILE(2,:) = ['swuc_wpi_AR1_02'];
DFILE(3,:) = ['swuc_wpi_AR1_03'];
DFILE(4,:) = ['swuc_wpi_AR1_04'];
DFILE(5,:) = ['swuc_wpi_AR1_05'];
DFILE(6,:) = ['swuc_wpi_AR1_06'];
DFILE(7,:) = ['swuc_wpi_AR1_07'];
DFILE(8,:) = ['swuc_wpi_AR1_08'];
DFILE(9,:) = ['swuc_wpi_AR1_09'];
DFILE(10,:) = ['swuc_wpi_AR1_10'];
DFILE(11,:) = ['swuc_wpi_AR1_11'];
DFILE(12,:) = ['swuc_wpi_AR1_12'];
DFILE(13,:) = ['swuc_wpi_AR1_13'];
DFILE(14,:) = ['swuc_wpi_AR1_14'];
DFILE(15,:) = ['swuc_wpi_AR1_15'];
DFILE(16,:) = ['swuc_wpi_AR1_16'];
DFILE(17,:) = ['swuc_wpi_AR1_17'];
DFILE(18,:) = ['swuc_wpi_AR1_18'];
DFILE(19,:) = ['swuc_wpi_AR1_19'];
DFILE(20,:) = ['swuc_wpi_AR1_20'];
DFILE(21,:) = ['swuc_wpi_AR1_21'];
DFILE(22,:) = ['swuc_wpi_AR1_22'];
DFILE(23,:) = ['swuc_wpi_AR1_23'];
DFILE(24,:) = ['swuc_wpi_AR1_24'];
DFILE(25,:) = ['swuc_wpi_AR1_25'];
DFILE(26,:) = ['swuc_wpi_AR1_26'];
DFILE(27,:) = ['swuc_wpi_AR1_27'];
DFILE(28,:) = ['swuc_wpi_AR1_28'];
DFILE(29,:) = ['swuc_wpi_AR1_29'];
DFILE(30,:) = ['swuc_wpi_AR1_30'];
DFILE(31,:) = ['swuc_wpi_AR1_31'];
DFILE(32,:) = ['swuc_wpi_AR1_32'];
DFILE(33,:) = ['swuc_wpi_AR1_33'];
DFILE(34,:) = ['swuc_wpi_AR1_34'];
DFILE(35,:) = ['swuc_wpi_AR1_35'];
DFILE(36,:) = ['swuc_wpi_AR1_36'];
DFILE(37,:) = ['swuc_wpi_AR1_37'];
DFILE(38,:) = ['swuc_wpi_AR1_38'];
DFILE(39,:) = ['swuc_wpi_AR1_39'];
DFILE(40,:) = ['swuc_wpi_AR1_40'];

% variables to be saved
% variables to be saved
varname(1,:) = ['SD']; % states
varname(2,:) = ['vD']; % transient volatility and signal noise process 
varname(3,:) = ['WD']; % covariance matrix for log volatility innovations
varname(4,:) = ['MD']; % standard deviation for measurement error
varname(5,:) = ['AD']; % autoregressive parameter for measurement error

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read data

% warren-pearson data 1749-1890 (1782-1784, 1788, 1792 are missing) 1914 = 100 
% Historical Statistics of the United States
% read data for 1797-1890
wp = xlsread('Cc113-124','B45:B138'); 
inf_wp = diff(log(wp));
T_wp = size(wp,1);
year_wp = 1797 + [0:1:T_wp-1];
% figure
% plot(year_wp(2:T_wp),inf_wp,'-b','Linewidth',2); hold on;
% axis([1795 2015 -inf inf])

% hanes data 1860-1990 (1942-46 are missing), 
% Historical Statistics of the United States
% For the period 1860-1941, 1914 = 100
hanes1 = xlsread('Cc125-137','B3:B84'); 
inf_h1 = diff(log(hanes1));
T_hanes1 = size(hanes1,1);
year_hanes1 = 1860 + [0:1:T_hanes1-1];
% plot(year_hanes1(2:T_hanes1),inf_h1,'-g','Linewidth',2); hold on;

% For the period 1947-1990, 1990 = 100
hanes2 = xlsread('Cc125-137','B86:B129'); 
inf_h2 = diff(log(hanes2));
T_hanes2 = size(hanes2,1);
year_hanes2 = 1947 + [0:1:T_hanes2-1];
% plot(year_hanes2(2:T_hanes2),inf_h2,'-g','Linewidth',2); hold on;

% bls data 1890-1997 (1982 = 100)
% Historical Statistics of the United States
bls_hs = xlsread('Cc66-83','B2:B109'); 
inf_bls_hs = diff(log(bls_hs));
T_bls_hs = size(bls_hs,1);
year_bls_hs = 1890 + [0:1:T_bls_hs-1];
%plot(year_bls_hs(2:T_bls_hs),inf_bls_hs,'-r','Linewidth',2); hold on;

% bls data from FRED (1982 = 100), monthly, NSA, 1/1913-6/2013
bls_m_fred = xlsread('PPI_ACO_FRED_July_16_2013','B13:B1218'); % monthly
Tm = size(bls_m_fred,1);
bls_fred = bls_m_fred(6:12:Tm); % annual, point sampled in June (maximizes correlation with annual BLS HSUS data)
inf_bls_fred = diff(log(bls_fred));
T_bls_fred = size(bls_fred,1);
year_bls_fred = 1913 + [0:1:T_bls_fred-1];
% plot(year_bls_fred(2:T_bls_fred),inf_bls_fred,'-m','Linewidth',2); hold off;
% legend('Warren-Pearson','Hanes 1891-1941','Hanes 1948-1990','BLS HSUS','BLS FRED')
% set(gca,'Fontsize',16)

% noisy inflation data Warren-Pearson 1798-1860,Hanes 1861-1941, BLS HSUS
% 1942-1947, Hanes 1948-1990
yn = [inf_wp(1:63); inf_h1; inf_bls_hs(52:57); inf_h2]; % noisy inflation data 
% clean inflation data
yc = inf_bls_fred(35:T_bls_fred-1); % BLS FRED 1948-2013
year = [1798:1:2013]';
T = size(year,1);
% figure
% plot(year(1:size(yn,1)),yn,'-g','Linewidth',2'); hold on;
% plot(year(151:T),yc,'-b','Linewidth',2'); hold off;
% legend('Noisy','Clean')
% set(gca,'Fontsize',16)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% partition sample
Y0 = yn(1:52); % training sample 1798-1849

year = [1850:1:2013]'; % years to be used for estimation
T1 = 98; % 1947
T2 = 141; %1990
T = size(year,1);

YS1 = yn(53:size(yn,1))'; % 1850-1990 to be used for estimation
YS2 = [zeros(1,T1) yc']; % 1948-2013 padded w/ zeros to mkae dating easier

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a weakly informative prior for initial states (ordered \pi, \mu)
SI = [mean(Y0)*ones(2,1); 0]; % prior mean on initial state
PI = [.15^2 0 0; 0 0.025^2 0; 0 0 .15^2]; % prior variance on initial state
R0 = var(Y0); % prior variance for SW transient innovations 
Q0 = R0/25; % prior variance for trend innovations 
% prior variance for log R0, log s0 (ballpark numbers)
ssr0 = 5;
ssq0 = 5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prior for W (inverse wishart w/ scale matrix W0 and degrees of freedom v0) (innovation variance for log volatilities)
v0 = 10;
sv0 = 0.2236*sqrt((v0+1)/v0); % stock and watson's calibrated value adjusted for time aggregation
TW0 = v0*(sv0^2)*eye(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prior for measurement-error variance \sigma_m 
vm0 = 2; % proper but weakly informative
sm0 = 0.5*sqrt(R0)*sqrt((vm0+1)/vm0); %mode is half of training sample standard deviation
dm0 = vm0*(sm0^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prior for measurement-error ar parameter \rho 
rho0 = 0; % prior mean on ar coefficients
Vrho0 = 0.45^2; % prior variance on ar coefficients
invVrho0 = inv(Vrho0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize gibbs arrays
SA = zeros(NG,3,T); % draws of the state vector
vA = zeros(2,T+1,NG); % stochastic volatilities, ordered (r,q) for SW+ parameterization
WA = zeros(NG,2,2); % standard error for log volatility innovations
SM = zeros(NG,1); % standard error for innovation to measurement error
AA = zeros(NG,1); % ar parameter for measurement error

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize stochastic volatilities and measurement error variance
vA(1,:,1) = R0*ones(T+1,1);
vA(2,:,1) = Q0*ones(T+1,1);
WA(1,:,:) = (sv0^2)*eye(2);
SM(1,:) = sm0;
AA(1,1) = rho0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% begin MCMC
%warning('off','all');
for file = 1:NF,
    file
    for iter = 2:NG,
        
       % states conditional on hyperparameters and data
       [S0,P0,P1] = kf_SWR(YS1,YS2,vA(2,:,iter-1),vA(1,:,iter-1),SM(iter-1),AA(iter-1),SI,PI,[T1,T2,T]);
       SA(iter,:,:) = GIBBS1_SWR(S0,P0,P1,AA(iter-1),T);
  
       % state innovations
       Saug = [SI squeeze(SA(iter,:,:))];
       f = ([1 -1 0; 0 1 0; 0 0 1]*Saug(:,2:T+1) - [0 0 0; 0 1 0; 0 0 AA(iter-1)]*Saug(:,1:T))'; 
       
       % stochastic volatilities
       vA(:,1,iter) = bivariate_svmh0(vA(:,2,iter-1),[R0;Q0],squeeze(WA(iter-1,:,:)),ssr0,ssq0);
       for t = 2:T, 
          vA(:,t,iter) = bivariate_svmh(vA(:,t+1,iter-1),vA(:,t-1,iter),squeeze(WA(iter-1,:,:)),f(t-1,1:2)',vA(:,t,iter-1)); 
       end
       vA(:,T+1,iter) = bivariate_svmhT(vA(:,T,iter),squeeze(WA(iter-1,:,:)),f(t-1,1:2)',vA(:,T+1,iter-1)); 
      
       % covariance matrix for log volatility innovations (inverse wishart draw)
       log_vol_inn = diff(log(vA(:,:,iter))');
       TW = TW0 + log_vol_inn'*log_vol_inn; % posterior scale matrix
       %PS = real(sqrtm(TW\eye(2,2))); %sqrtm is too slow
       [V,D] = eig(TW);
       TWsqrt = real(V*(D.^.5)*V');
       PS = TWsqrt\eye(2,2); 
       u = randn(2,v0+T);
       WA(iter,:,:) = (PS*u*u'*PS')\eye(2,2);
  
       % measurement error
       % innovation variance
       v = ig2(vm0,dm0,f(:,3)); 
       SM(iter,:) = v^.5;
       
       % autoregressive parameter 
       yy = squeeze(SA(iter,3,2:T))/SM(iter,:);
       XX = squeeze(SA(iter,3,1:T-1))/SM(iter,:);
       post_var = inv(invVrho0 + XX'*XX);
       post_mean = post_var*(invVrho0*rho0 + XX'*yy);
       AA(iter,:) = (post_mean + sqrtm(post_var)*randn(1,1))';
       
    end
    
    SD = SA(1:10:NG,:,:);
    vD = vA(:,:,1:10:NG);
    WD = WA(1:10:NG,:,:);
    MD = SM(1:10:NG,:);
    AD = AA(1:10:NG,:);
        
    save(DFILE(file,:),varname(1,:),varname(2,:),varname(3,:),varname(4,:),varname(5,:))

    % reinitialize gibbs arrays 
    SA(1,:) = SA(NG,:); 
    vA(:,:,1) = vA(:,:,NG); 
    WA(1,:,:) = WA(NG,:,:);
    SM(1,:) = SM(NG,:); 
    AA(1,:) = AA(NG,:);
   
end
exit


  